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ABSTRACT 


Many current satellites employ on-off thrusters to accomplish attitude control tasks 
which may include initial acquisition, rotational maneuvers, and on-orbit stabilization. This 
work shows that the use of pulse-width pulse-frequency (PWPF)-modulated thrusters 
provides several important advantages over conventional bang-bang thruster control 
methods, including less thruster activity and closer-to-linear actuation. The PWPF 
modulator is implemented in simulations using the Matrixx/Systembuild software package. 
Simulations assuming a rigid spacecraft are first performed to compare the performance of 
the PWPF-modulated thrust controller with that of conventional bang-bang and time- 
optimal bang-bang controllers. The discussion is then extended to the case of a spacecraft 
with structural flexibility, as is encountered quite often in three-axis stabilized vehicles with 
large fold-out solar arrays. Simulations for comparison of the controllers are performed 
using the flexible spacecraft dynamics model. The control loop design in the presence of 
flexibility and possible interaction with the PWPF modulator nonlinearity are addressed. 
Using a describing function model of the modulator, stability margin with respect to the 
structural mode limit cycle is predicted. Simulations are then conducted to verify the 
predicted stability margin. 
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L INTRODUCTION 


A. REASON FOR RESEARCH 

This work was performed in support of ongoing research in flexible spacecraft 
dynamics and control being conducted under the direction of Dr. Brij Agrawal. 
Justification is presented for implementing pulse-width pulse-frequency (PWPF) 
modulated thruster control as an alternative to conventional bang-bang control methods. 
The goal was to provide smoother control for improved pointing accuracy with less thruster 
activation and reduced excitation of the elastic modes of vibration. Simulation and analysis 
are presented to illustrate the improvement afforded by PWPF modulation and to address 
the stability issues which arise in the presence of nonlinear actuation. The research 
provides foundation for future thruster control experiments using the Spacecraft Dynamics 
Laboratory facilities. 

B. BACKGROUND 

1. Attitude Control System Overview 

A satellite attitude control system (ACS) is employed to provide vehicle 
acquisition, rotation maneuvers, and on-orbit stabilization. In acquisition, the control 
system acts to null initial angular velocity errors and bring the spacecraft to its operational 
configuration. This situation is often encountered immediately after spacecraft separation 
from the final launch vehicle stage, whereby the separation event has introduced initial 
angular rates. Rotational (or slewing) maneuvers are sometimes performed to reorient the 
spacecraft to meet mission objectives, or to prepare for translational propulsive maneuvers. 
On-orbit stabilization refers to the nominal mode of operation whereby the control system 
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typically acts in the regulator mode to maintain a particular pointing attitude with respect to 
the earth. 

An ACS must provide stability under the influence of various types of external 
disturbance torques. These may include the secular disturbance due to thrust misalignment 
during a translational propulsive maneuver, or cyclical environmental torques, e.g. from 
solar pressure or gravity gradient. Internal disturbances are applied by the motion of 
satellite components such as steerable antennas, sensor suites or sun-tracking solar arrays. 

2. Introduction to Thruster Attitude Control 

Most modem satellite attitude control systems incorporate several types of 
actuators to include both internal momentum exchange devices (e.g. momentum wheels, 
reaction wheels, control moment gyros) along with a complement of thrusters. Thrusters 
are used in situations when disturbance torques exceed the control authority of the 
momentum exchange devices, and are often capable of much faster reorientation 
maneuvers. PWPF-modulated thrusters have been used to maintain stability during 
translational maneuvers [Reference 1]. 

Momentum exchange type actuators can usually be treated as proportional 
devices. On the other hand, thrusters are typically on-off devices, only capable of 
providing a single fixed thrust level. Because this is inherently nonlinear actuation, thruster 
control systems cannot be analyzed in the same systematic and straightforward manner as a 
linear system. Response characteristics are not scalable as with a hnear system, but rather 
are highly dependent on input amplitude. By considering the system operation for the 
various attitude control modes (i.e., various levels of error signal amplitude), a general 
understanding of the control system behavior can be obtained. 

The main concern in design of a thruster control system is in performing the 
control operations within the required accuracy, with the least amount of thruster activity. 
Lower thruster activity directly translates into lower propellant consumption and less 
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thruster/ valve wear. Thus mission benefits such as longer mission life and/or more 
available payload weight for the same amount of fuel can be realized. 

3. The Complication of Flexibility 

Structural flexibility must be addressed in ACS design. In the presence of the 
flexibility, the possibility exists that natural frequencies of the structure can fall within the 
control bandwidth. This can lead to control-structure interaction problems which 
compromise the stability and pointing accuracy objectives (and possibly harm the 
structure). The dominant natural frequencies of vibration must be predicted, which often 
requires finite element modeling of the structure. The nonlinearity of the thruster-actuated 
control system can interact with the flexibility, resulting in self-sustained oscillation of the 
structure and back-and-forth thruster activity. 
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n. THRUSTER ATTITUDE CONTROL SYSTEM DESCRIPTION 


A. SIMPLE BANG-BANG CONTROL SYSTEM 

The block diagram representation of a simple satellite bang-bang control system with 
angular position and rate feedback is shown in Figure 1. The satellite, depicted in 
Figure 2, is modeled as a rigid rotational inertia which is controlled by a complement of 
thrusters which provide either a positive or negative torque of constant magnitude U or - U 
[adapted from Reference 2]. 


controller & thrusters plant (rigid satellite dynamics) 



Figure 1. Bang-Bang Attitude Control System 


Here 5 is the Laplace variable, I is the 
total spacecraft inertia, 9 is the controlled 
satellite pointing angle, and the factor 
]/Is2 is the familiar s-domain 
representation of the rigid body 
dynamics. The control torque of 
magnitude +/-U is commanded based on 
the sign of the error signal, e(t), which is 
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a combination of gained velocity and position error. For this system the controller is 
defined by 


\U sgn{e{t)) \i\e{t)\>a 

10 if |e(Ol < o. 


where 

^(0 = Qref - ^(0 - ad{t) 


( 1 ) 


and a is the deadband half-width 


A phase plane trajectory for the control system of Figure 1, ignoring the deadband for 
now, is shown in Figure 3. Since rate feedback introduces damping to the system, the 
phase plane trajectory shows that the time response decays toward the origin where both 
rate and position are zero. In the absence of rate feedback the system is marginally stable, 
with undamped oscillation occurring with a unique amplitude for each unique initial 
condition. 



Figure 3. Phase Trajectory for Damped Thruster ACS with No Deadband 
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Note that the switching line with no deadband passes through the origin of the phase 
space. Each period of constant thrust, from starting point A to switching point B, and from 
point B to switching point C, etc., forms a parabola in the phase space [Reference 3]. 

1. The Limit Cycle 

A limit cycle can be described as a self-sustained oscillation condition which 
arises in nonlinear systems. Viewed in terms of the phase plane trajectory, a limit cycle is 
a closed path which is approached from a starting condition either from inside the closed 
path (usually with the exception of the origin), or from the outside. The system trajectory 
will continue to trace the limit cycle path ad infinitum (given that it is a stable limit cycle). 
In thruster control, the "rigid body" limit cycle is due to the minimum on-time of the 
thruster (the "structural mode" limit cycle due to the interaction of the nonlinearity with a 
flexible structure is discussed in Chapter V). Consider the phase trajectory of Figure 4, 
the trajectory near the origin of the phase space for the bang-bang controller of Figure 1. 
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The deadband affords periods of coasting which are represented as horizontal 
lines in the phase plot. Each vertical segment (actually a small section of a parabola) 
corresponds to a thruster pulse. Due to the minimum on-time of the thrust, the minimum 
impulse and hence minimum change in rate, Ad^^, produces the ladder-step approach 
toward the origin. It is this minimum on-time or minimum impulse of the thrusters which 
produces the limit cycle seen near the origin. Note that a trajectory starting from inside the 
limit cycle boundary will also eventually end up in the limit cycle condition (assuming a 
nonzero rate). The frequency of the limit cycle is of great concern to the thruster control 
system designer because it dictates the fuel usage for the nominal mode of operation which 
is the predominant mode through the life of the satellite. However, the desire for a 
minimum frequency limit cycle is often in conflict with the pointing requirements since a 
tighter requirement on 0 for instance will result in less coast period and hence higher limit 
cycle frequency. Note that in Figure 3, after switching point C the absence of the 
deadband and associated coast periods cause "chatter" in the approach to the origin. This is 
wasteful of fuel, and upon arrival in the region of the origin, the limit cycle is of much 
greater frequency than in the system with the deadband. 

2. Time-Optimal Bang-Bang Control 

Time-optimal bang-bang control for a rest-to-rest maneuver is accomplished using 
the control law 


u(t) = -Usgnl 


le 


e 


^ ^cmd+ 2U 


( 2 ) 


where all terms are defined as in Figure 1. 

For any initial condition the phase plane trajectory consists of two parabolas with a 
single switching point, corresponding to full-on acceleration until the half-time of the 
maneuver, followed by full-on deceleration along the second parabolic trajectory straight to 
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the origin. This maneuver has been studied and modified in many reports, since it is the 
most time-efficient approach. 

B. PULSE-WIDTH PULSE-FREQUENCY MODULATION 

Most thrusters used for attitude control are designed for on-off operation. Pulse-width 
pulse-frequency (PWPF) modulation provides a means of producing a variable average 
thrust using a simple on-off thruster, hence an approximation to proportional control is 
achieved. 

A PWPF modulator is depicted in Figure 5. The average output of the pulse train is 
approximately proportional to the demanded torque input. Other forms of pulse 
modulation, for instance pulse width (or pulse frequency) or derived rate modulation are 
used in other applications which are not considered here. 



The modulator includes a relay with dead-zone and hysteresis, or Schmidt trigger, and 
an integrator (filter). A limiter after the filter was included in the digital implementation. 
The feedback from the output of the relay is subtracted from the input signal A 2 . The 
formation of pulses can be explained with the aid of Figure 6, which depicts the output of 
the discretized modulator and the filter used in this study, for a constant input. The filter 
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integrates the error signal (A 2 -A 3 ). Before the formation of a pulse the feedback from the 
relay is zero. So the magnitude of A 2 determines the rate of integration, and thus the rate 
of growth of the output of the filter. When the filter output reaches the threshold Con of the 
Schmidt trigger, the relay switches on and turns on the thruster. The feedback signal then 
resets the integrator at a rate proportional to (Um-Ai) until the filter output falls below eojf 
at which time the thruster signal is turned off and the feedback signal is again zeroed out. 
Note that the appearance of the pulses as being other than square is due to the integration 
time step of 10 milliseconds. The first pulse is commanded on at .06 seconds and 
commanded off at .08 seconds. 



0 .05 .1 .15 .2 .25 3 .4 .45 .5 

thne^ac. 


Figure 6. PWPF Modulator Output for Constant Input 
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Reference 1 notes the modulator tends to contribute control loop gain reduction at low 

frequency and phase lag at higher frequency, and presents design parameters based on 

tradeoffs between phase lag and internal deadband size: 

^.= 4.5 
r =0.12976 

m 

e,„=0.45 (3) 

The modulator was discretized for digital simulation. Using the signal labels from 
Figure 5, the filter output can be written 

Considering the signal at two consecutive time steps n and (n+1), and replacing the Laplace 
variable with its time equivalent (a discrete-time first derivative), Equation (4) can be 
reformulated as 

(Aj - A^)K = Ai+ sA^T 

Solving for the new value of the filter output, 

where dt is the sampling interval of the digital simulation to represent the modulator, which 
is chosen to produce a minimum on-time or minimum pulse width of dt. The digital 
implementation of the modulator is included in Appendix A. Simulation of the PWPF 
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modulator was performed using a sampling interval of 10 milliseconds. Figure 7 
illustrates operation of the modulator for two sinusoidal inputs and a ramp. 



Figure 7. PWPF Modulator Output for Various Inputs 



The modulation factor or duty cycle is defined as the average output of the 
modulator. Pulse trains for various constant inputs were analyzed to verify the expected 
linear relationship between input and duty cycle. The pulse trains are included in 
Appendix A. Once the pulse width and frequency are determined, the duty cycle is 
simply determined by: 

duty cycle = pulse width x pulse frequency (7) 

The output characteristics are tabulated in Table 1 and shown graphically in Figure 8. 
The gains computed in the table are the average values of the duty cycle for a given input. 
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TABLE 1. MODULATOR RESPONSE TO CONSTANT INPUTS 



Figure 8. Output Pulse Width, Pulse Frequency, and Duty Cycle 
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The modulation factor versus input graph of Figure 8 verifies the nearly proportional 
operation. Note that the modulator has an internal deadband, whereby a pulse will not be 
triggered for an input signal lower than about 0.15. It is important to remember that even 
though the output relationship is fairly linear, the relay constitutes a hard nonlinearity and 
hence nonlinear analysis must be employed to study the details of control system 
performance. 


C. INTRODUCTION TO SIMULATIONS 

Digital simulations were used for verifying basic control system operation, comparing 
the performance of different types of thruster control, and exploring system response with 
structural flexibility and variation of various parameters. In most cases simplified block 
diagrams are included in the text portion of the thesis for clarity and ease of interpretation, 
with corresponding Matrixx/Systembuild block diagrams included in Appendix B. All 
simulations were performed for the single axis control case, neglecting cross-axis 
momentum coupling effects. Hence the rigid body equation of motion is simply: 

Id = r, (8) 


from which the Laplace transform yields the s-domain transfer function representation: 


0(5) ^ 1 

TM Is^ 


(9) 


where 6 is the controlled pointing angle, I is the total moment of inertia, and Tc the 
control torque. For the case of thruster control Tc is approximated as a square pulse. 

Unless otherwise noted, all block diagram elements were discretized with a sampling 
interval of 10 milliseconds. This enforced the minimum on-time of the thruster, while 
providing a sample rate fast enough to sufficiently capture the dynamic response. 
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All simulations were run using the parameters of the flexible spacecraft simulator 
(FSS) experimental system (but assumed a rigid body for the simulations of Chapter III). 
The values used for total moment of inertia and thruster control torque were 11.4 kg-m^ 
and 0.35 N-m, respectively. The FSS was designed as a scale version of the pitch axis of 
a typical geosynchironous momentum-biased satellite, so that results representative of a 
realistic satellite control system could be expected. Indeed these values scale nicely to the 
pitch axis inertia and control torque values (440 kg-m^ and 10 N-m) for the operational 
satellite of Reference 1. 

A pointing accuracy goal of +/- .1 degrees was used throughout the study as a basis 
for comparison. This is a realistic (if not lenient) requirement for a communications 
satellite at geosynchronous altitude. 
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in. COMPARISON OF THRUSTER CONTROL METHODS FOR RIGID 

SPACECRAFT 


Simulations were run to demonstrate the performance of the pulse modulated control 
system and to compare it with the performance of bang-bang control systems. A block 
diagram control system representation for rigid body simulations is shown in Figure 9. 
Implementation of the pulse modulator involved replacing the elements from Figure 5 with 
a single Usercode Block in Systembuild which calls the FORTRAN routine for the discrete 
modulator in each integration step. For the time-optimal bang-bang control system the 


CONTROLLERS 



Figure 9. Attitude Control Systems for Rigid Body Simulations 


15 





















feedback gain and forward loop gain are not applicable. A linear (non-thruster) controller 
with position and rate feedback would bypass the controllers and the thruster block. The 
PWPF modulator incorporates an additional gain K which will be considered separately 
from the gain ki. The disturbance torque is included in some simulations. The error 
signal, e(t ), includes the error between the commanded angular position and the actual 
angular position, and the error between the commanded and actual angular rates where the 
commanded rate is always implicitly zero. 

The linear controller was used to size gains for a particular desired performance. The 
same gains were then used with the nonlinear controllers and for the simulations with 
flexible dynamics, in order to form a basis for comparison. The closed loop transfer 
function of the system for the case of linear actuation is: 




kjj 


2 k,X h 
5^ -|--^5-hy 


( 10 ) 


The equation is of the same form as the general equation for a viscously damped 
single-degree-of-freedom system. 


2- 2 

s + 2 ^ 0 }„s + co„ 

where and cOn are the damping ratio and natural frequency, respectively. Equating 
coefficients of (10) and (11) it follows that 

and C = (12) 

For the linear design a damping factor of .84 was chosen for minimal overshoot. The 
value of Q)n can be thought of as the control bandwidth, and was selected as .06 Hz in 
order to minimize interaction with the first pole and zero of the flexible dynamics in later 
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simulations (.12Hz and .15Hz, respectively). Based on these values Equations (12) yield 
the gain values 

ki = 1.4 N-m/rad 
kiX= 6.7 N-m-s/rad 

The response to a ten degree step input with linear actuation is shown below in 
Figure 10. Graphs are shown for angular position in degrees, angular rate in degrees per 
second, the gained error signal kje(t), and the phase plane trajectory. A 2% settling time of 
13.6 seconds is achieved, with less than one percent overshoot. 



Figure 10. Step Response for Linear Actuation with Position and Rate 

Feedback 
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A. THRUSTER CONTROL STEP RESPONSE (SLEW MANEUVER) 

For large values of the error signal the system response with thrust-modulation is 
similar to the simple bang-bang system because the PWPF modulator tends to operate 
above its linear range. This is demonstrated by comparing Figures 11 and 12 below, 
where a 10° slew maneuver has been commanded from rest for bang-bang and modulated 
thruster control. Both controllers complete the maneuver in about 24 seconds, slower than 
the linear case but with essentially no overshoot as was desired. 



Figure 11. Slew Maneuver with Bang-Bang Controller 


The control torque pulses change in frequency for both controllers. This can be 
understood by referring back to Figure 4, which depicted the "ladder-step" phase plane 
approach to the final commanded state. The coast period (horizontal line in phase plane) 
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between each pulse increases as the trajectoiy approaches 0 = 0 , hence the pulse frequency 
decreases. The trajectory rides the boundary of the deadband from just after one second 
until the completion of the maneuver. 



Figure 12. Slew Maneuver with Thrust-Modulated Controller, K=39 


Note that the modulated maneuver is slightly smoother in the phase plane. The 
response is closer to the characteristically smooth phase trajectory of the linear system 
(Figure 10) because the PWPF modulator is closer to approximating linear actuation. The 
smoothness of the trajectory is an important factor in the presence of flexibility. An abrupt 
change in velocity curve from positive slope to negative slope indicates an abrupt change in 
acceleration and hence jerk as well. The faster the acceleration is reversed, the more jerk is 
applied, which for a flexible structure results in more stored flexural energy. Whereas the 
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bang-bang controller only operates according to the sign of the error signal, the PWPF 
modulator responds to both the magnitude and the sign of the error. The difference in 
phase trajectory smoothness is shown more dramatically in simulations where the feedback 
gains and kj are all set to unity (Figures 13-15). The resulting damping ratio of .15 (for 
the linear controller) produces excessive overshoot but also a faster rise time. Both of these 
effects are seen to some extent with the nonlinear controllers, however the modulated 
system provides a smoother transition between the reversal of acceleration. 



Figure 13. Underdamped Slew Maneuver with Linear Controller 
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Figure 14. Underdamped Slew Maneuver with Simple Bang-Bang Control 



Figure 15. Underdamped Slew Maneuver with PWPF Modulator, K=30 
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Note from Figure 15 that the pulse-modulated response would closer match the 
linear simulation if the modulator were operating in its linear region over the maneuver. As 
the time history of gained error signal shows, the modulator input is above the modulator 
saturation level of 1.05 (see Figure 8) until about 2.5 seconds, causing continuous-on 
operation. 

For the linear control loop, application of the final value theorem indicates that the 
steady state error to a step input is zero. The PWPF modulator on the other hand contains 
an internal dead zone which results in a steady state error (or defines the amplitude of the 
rigid body limit cycle), such as can be seen in Figure 15. The small steady state error is 
achieved only by using a high value of modulator input gain K (=39 for the previous runs). 
The need for a substantial value of input gain can be understood by considering the 
magnitude of the signal prior to the modulator input gain. The initial error signal is due to 
the commanded input only, since the spacecraft is initially at rest: 

e(0)x^= X 1.0 = .174 

which is the maximum value of this signal over the entire maneuver, already quite small 
considering the modulator dead zone of approximately +/- 0.15. Because the thruster is 
full-on for the first 2.5 seconds, the angular rate increases quickly to nullify the error 
signal: 

e(~2.5sec)xA:, = (10°x 1 - 6°xl - 4.2°/sec xl) x ^= -.003 

which is much too low to trigger the modulator. So a dramatic gain factor is needed to 
keep the modulator input signal above its internal deadband. It was found, however, that 
the value of modulator input gain could be substantially lower, say K =10, while still 
achieving the maneuver (Figure 16). 
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Figure 16. Underdamped Slew Maneuver with PWPF Modulator, K=10 


The rise time is not affected appreciably by the lower input gain, since the modulator is 
still initially saturated, producing a period of constant thrust. The lower value of K 
however results in more steady state error because the deadband is reached sooner. As 
with a linear system, an increase in the loop gain (via the modulator input gain) decreases 
the steady state error in the presence of a disturbance torque; this is discussed further in 
Chapter IV. 

Note that limit cycle behavior was not seen after maneuver completion in the previous 
simulations. This is because the final thruster pulse brings the velocity to zero. With non¬ 
zero initial conditions the limit cycle is seen at the completion of the maneuver (Figure 
17). 
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Figure 17. Thrust-Modulated Slew Maneuver with Non-Zero Initial 

Conditions 

A ten degree slew maneuver performed using the time-optimal bang-bang controller is 
shown in Figure 18. Although the maneuver is faster than the other types of thruster 
control, the acceleration is by definition of the control law reversed instantaneously at the 
half-maneuver time. The effects of this with a flexible structure can be seen in the results 
presented in Reference 4, which show the maximum strain in the structure occurring just 
after the half-maneuver point. The stored energy then produces a whip-like motion of the 
flexible appendage which contributes to overshoot at the end of the maneuver. In rigid 
body simulations without dead zone, the sampling interval alone constituted enough of a 
delay to cause switching instability at the end of the maneuver. Hence it was necessary to 
include dead bands on position error, {6^^ - 6^,), and rate , in order to avoid thruster 

chatter at the sampling rate of 100 Hz. Even with the dead zone. Figure 18 shows a slight 
imperfection in the half-time switching. Note that the maximum velocity reached 
(approximately 4.5 deg/sec) is over twice the velocity achieved using the other controllers. 
In many systems rate limiting is necessary, which of course diminishes the intended 
advantage of time-optimal switching. 
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Figure 18. Time-Optimal Bang-Bang Maneuver 


Sensitivity of the time-optimal control system to parameter uncertainty was studied. 
Note that the switching law (Equation (2)) contains the thrust value and the total moment of 
inertia. Since these parameters are generally not known exactly, the values used in the 
switching law are only a best estimate. Though a thruster may be accurately calibrated, the 
actual torque produced is subject to mounting misalignment tolerances and plume 
impingement effects. Also, the moment of inertia of the satellite changes over the orbit 
lifetime with the depletion of propellant. Uncertainty manifests as overshoot (or 
undershoot) and subsequent unintended thruster firing (see Figures 19 and 20). Note 
the rigid body limit cycle which is established after completion of the maneuver. 
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Figure 19. Time-Optimal Bang-Bang, 20% Overestimate of Moment of Inertia 


Figure 20. Time-Optimal Bang-Bang, 20% Underestimate of Moment of Inertia 
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With the PWPF-modulator uncertainty does not play as critical a role since an estimate 
of the inertia and the control torque are not part of the control law formulation. 

B. NOMINAL CONTROL MODE: RIGID BODY LIMIT CYCLE 
PERFORMANCE 

The rigid body limit cycle performance was determined for each controller by 
simulating the response to small initial conditions. Limit cycles of the same size in the 
phase plane correspond to equivalent limit cycle frequencies. With the goal of keeping 
within +/- 0.1° pointing accuracy, the deadband of the bang-bang controller and the input 
gain for the modulated controller were adjusted accordingly. In the case of the simple on- 
off controller with deadband, the limit cycle behavior is shown in Figure 21. From the 
position graph the limit cycle frequency is determined to be about .026 Hz. 



Figure 21. Rigid-Body Limit Cycle with Bang-Bang Controller 
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A limit cycle of similar frequency is achieved with the PWPF modulator (Figure 22) 
with an input gain K of 39. Treatment of the rigid body limit cycle in more general 
mathematical terms for any thruster system can be found in Reference 5, which also notes 
that with position and rate feedback, the limit cycle is stable. 



Figure 22. Rigid-Body Limit Cycle with Thrust-Modulated Controller 


The same limit cycle is established for both controllers because minimum-impulse 
operation is established, i.e. each vertical portion of the limit cycle corresponds to the 
minimum .01 second thruster pulse (with .02 deg/sec). This is unexpected since 

early literature indicates a conventional bang-bang controller would not converge to 
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minimum impulse operation, but would rather form the higher frequency limit cycle shown 
in Figure 23 [adapted from Reference 6]. 



Figure 23. Expected Limit Cycles for Conventional (Bang-Bang) and 
Minimum Impulse (PWPF) Systems 

Though time-optimal switching would probably not be used for nominal attitude 
control, a limit cycle of comparable size can be produced by introducing plant inertia 
uncertainty (here -20% was used as in Figure 20) . Figure 24 shows that minimum- 
impulse operation is not established and the resulting frequency is determined to be about 
. 18 Hz, substantially higher than the other controllers. 
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Figure 24. Rigid-Body Limit Cycle with Time-Optimal Switching 


C. DISTURBANCE REJECTION 

The response of the various control systems to a constant disturbance torque was 
studied. A constant torque of .015 N-m was applied in most cases. This value was 
originally used in Reference 7 to model the torque applied due to cabling wind-up in the 
experimental FSS setup. The value scales nicely to the thruster misalignment disturbance 
torque countered by the controller in Reference 1. The assumption of a constant 
disturbance torque is adequate for the case of disturbance due to thruster misalignment 
during a translational maneuver. The response to stochastic disturbances such as sensor 
noise are not considered. 
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A typical phase portrait for the thruster systems [Ref. 5, Ref. 6] converges to a 
disturbed limit cycle of the shape shown in Figure 25. The time history of Figure 26 
shows that as the deadband is reached due to the disturbance, a thrust pulse is commanded 
to counter the disturbance. The plot of 6 shows that the satellite periodically "bounces" off 
the deadband. These plots were generated using a more powerful thruster and are therefore 
intended for illustrative purposes only. 



Figure 25. Phase Plane Trajectory with Convergence to Disturbed 

Limit Cycle 



Figure 26. Typical Response to Constant Disturbance Torque 


Using the system parameters considered thus far, the angular position responses for 
the bang-bang and the PWPF-modulated systems (with K=39) are shown in Figures 27 
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and 28. Both results show thruster activity at 4.2 Hz, where the disturbance torque keeps 
the satellite at the edge of the deadband. 



Figure 28. Response of Pulse-Modulated Control System to Constant 


Disturbance Torque 

Note that the modulated controller does not provide the required pointing accuracy of 
.1 degrees. With a linear controller the steady state error under the influence of a constant 
disturbance torque is inversely proportional to the controller gain. With the PWPF 


32 













modulator the gain on the error signal can similarly be increased in order to lower the 
steady state error. This is demonstrated in flexible body simulations in Chapter IV. The 
flexible body stability problems which can be encountered due to the thruster activity and 
the increased gain are further discussed in Chapters IV and V. 

A closer look at the gained error signal and the resulting pulses (Figure 29) for the 
PWPF-modulated simulation shows that when the error signal reaches the internal 
deadband of the modulator (point A) a single pulse is fired which decreases the error signal 
until the end of the pulse (point B). The error signal then grows under the influence of the 
disturbance until the deadband is reached again, triggering another thruster firing (point C). 
The pulse width of ten milliseconds is consistent with the response to the input of -.15 as 
given in Table 1, however the 8 Hz pulse rate cannot be realized since the input signal is 
not constant, but rather drops with each pulse. 



Figure 29. Error Signal During Constant Disturbance Torque 
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IV. FLEXIBLE SPACECRAFT DYNAMICS AND SIMULATIONS 


The simulations in Chapter HI demonstrated the advantages of PWPF modulation for 
the rigid body case. The flexible body dynamics are now discussed and comparative 
simulations presented. 

A. FLEXIBLE SPACECRAFT MODELS 

It is often necessary to include several vibrational modes of a flexible structure in a 
spacecraft dynamical model in order to adequately design the control system and to foresee 
possible modes of interaction between the control system and the structure. It is 
instructional to first examine the model of a spacecraft with one dominant mode of vibration 
about the controlled axis. Consider the spacecraft with a flexible solar array pictured in 
Figure 30 [adapted from Reference 8]. 



Figure 30: Spacecraft With Flexibility 
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Assuming a single undamped dominant torsional mode of vibration, the transfer 
function from control input to angular position output can be represented as; 


G(5) = 


d,is) ^ s^ + (0„^ 
T(s) 


(13) 



and 


a=Q).J\ + ^ 


(14) 


where s is the Laplace transform variable, h and Is are respectively the rigid body and 
solar array inertias about the pictured axis, 0 i> is the controlled rotation angle of the rigid 
central body,7’ is the control torque, and K is the torsional stiffness. The term (0„ in the 
numerator is the cantilever frequency (fixed base frequency) of the flexible appendage, and 
Qn in the denominator is the "system mode" or "free-free" mode frequency of vibration, 
whereby the rigid hub is no longer fixed but vibrates together with the flexible appendage. 
The 52 temi in the denominator is attributable to the rigid body mode of the spacecraft (zero 
frequency). An important observation to be made is that the complex zeros of the plant 
transfer function correspond to the fixed-base vibrational frequency of the flexible 
structure. So the transfer function implies that while the flexible appendage may be 
vibrating at its natural frequency, a control torque is supplied to keep the central body from 
rotating [Ref 9]. For the more general case of N modes of vibration, the transfer function 
representation of the dynamics can be represented as 


G(s) = 


0 ( 5 ) 1 A [//m,^ + (2CM)5 + l] 

ns) /54=|u'/^i.'+(2C,M)5 + l] 


(15) 


where I is the entire spacecraft inertia, (Oi is the ith cantilever mode natural frequency, and 
Qi is the ith system mode natural frequency. 
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The dynamics of Equation (13) can be derived easily because of the simplification that 
a single vibrational mode and a known stiffness K could represent the solar array 
flexibility. A higher fidelity dynamical model of a flexible spacecraft can be developed via 
the hybrid-coordinate representation of the system [Ref. 10], which includes both physical 
coordinates and modal coordinates. A finite element analysis is performed to determine the 
cantilever frequencies and mode shapes of the structure. For the general case of three 
rotational degrees of freedom, the linearized equations of motion are [Ref. 1], for the rigid 
body: 

ie + Hd + gq = T^ + T^ (16) 

and for the flexible appendage: 

Q + 2.C (Op + Q DO = 0 (17) 

where 1 is the entire spacecraft inertia matrix, H is the momentum coupling matrix, D is the 
rigid-elastic coupling matrix of the flexible appendage, 6 is the attitude angle vector, Tc 
and Td are the control and disturbance torque vectors, q is the cantilever modal coordinate 
vector, andtu^ is a diagonal matrix of squared cantilever modal frequencies. The inclusion 
of the damping term of Equation (17) is based on the assumption of modal damping, 
whereby the damping matrix has been uncoupled and a single value for the modal damping 
ratio, is assumed for all modes. 

The flexible dynamics model used in this study was derived [Ref. 7, Ref. 11] for the 
experimental flexible spacecraft system (ESS) using the hybrid-coordinate formulation. 
The system, pictured in Figure 31, free-floats on air pads and vibrates in planar motion, 
with translational motion restricted by the air bearing. Hence the momentum coupling term 
H (Equation 16) does not apply and I, 6, Tc and 7^ become scalars. Equations (16) and 
(17) then become [Ref. 7]: 
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(15) 


i=l 

q.+2C.Q}.q. + CO,^g. + D.e = 0 

with the rigid-elastic coupling matrix having reduced to a vector with each element based on 
the modal vector components in the rigid body frame (<})■ and ) and defined for each 

mode i of the flexible appendage as; 

A = - yF<Pi )dm ( 16 ) 

Each Di is calculated using a FORTRAN subroutine which operates on the cantilever modal 
frequencies and mode shapes from the finite element analysis output; the integral of 
Equation (16) becomes a summation over the number of discrete subbodies into which the 
flexible body has been divided, where Xf and yp- locate each subbody in the rigid center 
body frame. 


THRUSTER SYSTEM 
AIR TANK 

MOMENTUM WHEEL 
ASSEMBLY 

MAIN BODY 


Figure 31. Flexible Spacecraft Simulator 

In discretizing the system via the finite element method, the number of modes was 
truncated at six to obtain a compromise between reasonable model accuracy and 
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computational feasibility. The model was placed into state space in preparation for digital 
simulation using the Matrixx/Systembuild software package. The state space representation 

of the system equations is: 


x = Ax+ Bu 
y = Cx + Du 

where the state vector, x, is defined as: 

x = [0 . Qf, B 4,. 


( 20 ) 


( 21 ) 


The output y is the vector of the current states, hence C is a 14x14 identity matrix and it is 
assumed that the feedback values of angular position and rate are measured exactly. The 
system matrix form used in Matrixx/Systembuild is 


S = 


'A 

B' 

C 

D 


( 22 ) 


The system matrix was discretized for simulations with dt= .01 second. The complete 
mathematical descriptions of the state matrix A and the input matrix B are given in 
Reference 11 (p. 20). The direct transmission matrix D is set to zero. 

During initial simulations with the flexible dynamics, the natural frequencies of 
oscillation were found to vary somewhat from the values presented in Reference 7 (pages 
21 and 110). In order to find the actual frequencies present in the discretized model, a bode 
plot of the system was obtained (Figure 32). The DBODE command in Matrixx was 
used and the normalized frequencies were divided by the sampling period to obtain rad/sec. 
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Figure 32. Bode Plot of Discretized Plant Model 


Based on the bode plot the bode plot the cantilever and system frequencies were 
determined and are given in Table 2. 


TABLE 2. FSS DISCRETE MODEL CANTILEVER AND SYSTEM 

FREQUENCIES 


Mode 

Cantilever Frequency 
(Hz) (rad/sec) 

System Frequency 
(Hz) (rad/sec) 

1 

0.121 

0.76 

0.153 

0.96 

2 

0.330 

2.07 

0.410 

2.58 

3 

2.660 

16.71 

2.730 

17.15 

4 

3.550 

22.31 

3.630 

22.81 

5 

6.140 

38.58 

6.210 

39.02 

6 

16.610 

104.36 

16.900 

106.19 
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Reference 7 determined the first cantilever natural frequency to be .139 Hz, and the 
first system frequency .175 Hz. The discrepancy in the frequencies could partly be due to 
the method of discretization used. Consider the following continuous representation of a 
second order notch filter [Ref. 1]: 




(23) 


Figure 33 shows the continuous frequency response plotted with the discrete-time 
response using two different discretization methods (with w— .96). Note the shift in zero 
frequency introduced by the zero order hold method, whereas the first order hold method 
provides much better results. The same shift may have occurred when the FSS dynamics 
model was discretized. (This illustration underscores the difficulty in designing a notch 
filter for lightly damped poles since only a small inaccuracy in the knowledge of the true 
pole location would result in poor filter performance.) 



(a) (b) 

Figure 33. Continuos to Discrete-Time Conversion of Notch Filter using 
(a) zero order hold method, (h) first order hold method. 













































The state space representation of the system contains the information to form transfer 
functions between the torque input and each of the fourteen outputs. The transfer function 
for the central body rotation angle output was given in Equation (15), and can be modified 
into the form 


rrct - ^ - J-nf iV "in [£!±2C6>£+^ 


Using the values of the cantilever and system natural frequencies from Table 2, and 
an assumed modal damping ratio of .4 percent, the transfer function for the first three 
modes becomes: 


23 * (s^ +■ 0065+. 58)(5^ +.0175 + 4.41)(5^ +. 134^ + 279) 
(sA.0085+.922)(5V.02U + 6.76)(5^ +1385 + 296) 


(25) 


The poles and zeros of the complete six-mode transfer function were computed using a 
Matlab routine which is included in Appendix C and are shown in Figure 34. 



Real Axis 


Figure 34: Poles and Zeros of Transfer Function from T to 6 
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Note that the underdamped modes appear as alternating poles and zeros near the 
imaginary axis. Reference 9 shows that with position and rate feedback or lead 
compensation, stability of the modes is achieved, that is the root locus from each open loop 
pole to its adjacent zero will lie in the left half plane. With angular position and rate 
feedback the loop transfer function is multiplied by the factor ki{ ts+1) . The resulting root 
locus plot near the origin is shown in Figure 35, for z = A.ll. An underlying assumption 
is that the sensor and actuator are colocated at the "rigid" central hub of the spacecraft. 
With flexibility between the sensor and actuator, lead compensation can actually cause 
structural mode instability (see Reference 12, Appendix A.4). 



Figure 35. Root Locus vs. Forward Loop Gain With Position 

and Rate Feedback 
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If the controller gain were chosen to locate the closed loop poles near points A for 
maximum damping of the flexible modes, the resulting rigid body control would be 
considerably underdamped. More realistically, the closed loop poles would probably be 
placed with the rigid body pole closer to point B, for nearly-critical damping of the center 
body. Hence the flexible poles would be dictated by the rigid body controller, and would 
be much nearer the imaginary axis. Although the stable root locus plot is guaranteed with 
the linear controller, nonlinearity can cause closed loop pole locations to change as a 
function of time. Nonlinear analysis is discussed in Chapter V. 
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B. COMPARISON OF THRUSTER CONTROL METHODS FOR 

FLEXIBLE SPACECRAFT 

Comparisons of bang-bang and pulse-modulated thruster control in the presence of 
structural flexibility are presented in this section. The first two of the six cantilever modal 
coordinates are plotted to give some feel for the response of the flexible appendage with 
respect to the central body. The simulations were run at the same sampling rate (100 Hz) 
and with the same basic block diagrams as the previous rigid body simulations, only the 
flexible body dynamics were substituted for the rigid. 

1. Step Response 

The ten-degree step response is shown in Figures 36-41 for the PWPF- 
modulator, bang-bang, and time-optimal bang-bang controllers. The time responses of the 
first two modal coordinates are included to provide some indication of the arm excitation. 
The PWPF-modulator input gain of 30 (Figure 36) was reduced to 10 for the run of 
Figure 37. Note the reduced thruster firing with little loss in accuracy, as well as the 
rapid decay in the modal coordinate amplitudes. The oscillation in pointing angle is due to 
the input gain to the modulator, which increases the loop gain. Hence the original linear 
design with low controller bandwidth to avoid interaction with the first mode is not 
realized. It is possible that further "optimization" of the modulator could help improve this 
aspect of the results. 

For the first bang-bang run (Figure 38), the dead bands were left as were 
required to provide the minimum-impulse limit cycle from the rigid body simulations (with 
+/-.1° excursion in 6 ). The resulting thruster actuation is unsatisfactory, although tight 

control of the center body is noted. The forward deadband was opened up to +/- .2° to 
produce the result of Figure 39. The increased deadband reduces the "density" of 
thruster firings, however the oscillations and firings continue until almost 200 seconds. 
Similar results were obtained using the time-optimal controller. The first run (Figure 40) 
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shows a much faster maneuver (as is expected), however excessive switching occurs at the 
termination. Again opening up the deadband (Figure 41) provides virtually no benefit in 
terms of reducing thruster firings, and results in undesirable oscillations in pointing angle. 
The phase plane plot is shown in this run to point out the limit-cycle-like behavior, 
however further investigation showed thruster firing ends after about 200 seconds. The 
delayed overshoot occurring just after six seconds is probably due to the "whipping" 
motion of the flexible appendage as it releases the energy stored during the maneuver. The 
amplitude of cantilever modal coordinate qj reaches .26, much higher than for the thrust- 
modulated maneuver where the amplitude is .085. The rate of decay is also more gradual 
than the modulator case. 



t(ma(aae) 


Figure 36. PWPF-Modulated Step Response, K=30, Flexible Dynamics 
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Figure 37. PWPF-Modulated Step Response, K=10, Flexible Dynamics 



Figure 38. Simple Bang-Bang Step Response, Flexible Dynamics 
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Figure 39. Simple Bang-Bang Response with Increased Deadband 



Figure 40. Time-Optimal Bang-Bang Step Response, Flexible Dynamics 
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Figure 41. Time-Optimal Bang-Bang Response with Increased Deadband 


Though these simulations were intended to point out the advantages of PWPF- 
modulated thruster control, it is important to remember that variations of each control law 
have been devised to improve the performance in the presence of flexibility [e.g. Agrawal 
& Bang, Dodds]. The basic idea in slew maneuvers is to minimize the post-maneuver 
energy in the flexible structure, however this desire tends to conflict with the fact that a fast 
slew is also desired. [Hailey] demonstrated a slewing method which developed a more 
gradual angular acceleration using torque-shaping and resulted in less excitation of the 
structure. This method and others would show similar improvement for any of the 
controllers. 
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2. Small Initial Condition Response 

Figures 42-45 show the responses of each of the control systems to small 
initial conditions, where phase plane trajectories are included rather than the modal 
responses. Thruster activity in the first PWPF case (Figure 42) was reduced further by 
implementing an external deadband of +/- .1° before the modulator (Figure 43). The 
tradeoff is seen to be slightly less accuracy, but still well within the goal of +/- .1°. The 
simple bang-bang response (Figure 44) is comparable to the modulated response. The 
time-optimal bang-bang, which implements the dead bands directly in the feedback loop, is 
seen to bounce back and forth within the +/- .1° deadband. 



Figure 42. PWPF-Modulated Response to Small Initial Conditions, 

Flexible Dynamics 


49 

















me«liiM«r output 



Figure 44. Simple Bang-Bang Initial Condition Response 
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Figure 45. Time-Optimal Bang-Bang Initial Condition Response 
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3. Disturbance Rejection with PWPF-Modulated Controller 

The behavior of the modulated control system during a constant disturbance is 
illustrated in Figures 46-50. The disturbance of .015 N-m was applied for 20 seconds 
(except for Figure 46), starting at one second. Figure 46 shows that with ^=30, the 
required pointing accuracy appears to be achieved. After the disturbance, a rigid body limit 
cycle of frequency .05 Hz appears to have been established. The longer run of Figure 47 
(70 second disturbance) clearly shows that the required pointing accuracy is maintained. 
The results of Figures 48-50 show that as the input gain is increased, the steady-state 
error is reduced, however some degree of modulator instability is shown for the case of K= 
200 (Figure 50). After the disturbance, the over-excited modulator still manages to drive 
the system toward zero pointing error. 



Figure 46. Modulated Controller Response to Disturbance Torque, K=30 
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ire 47. Modulated Controller Response to Disturbance Torque, K=30 



Figure 48. Modulated Controller Response to Disturbance Torque, K=50 
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V. NONLINEAR CONTROL LOOP STABILITY ANALYSIS IN THE 
PRESENCE OF STRUCTURAL FLEXIBILITY 


The stability of the control system is a function of all elements of the control loop. The 
various elements other than constant gain add phase lead or lag (along with gain 
contribution except in the case of pure time delay), and hence contribute to the overall gain 
margin and phase margin. As was previously stated, for a control loop with linear 
actuation and position and rate feedback, the stability of each mode is assured. Though the 
control loop may induce oscillations in the structure, the oscillations will experience 
exponential decay. The discontinuous nature of the PWPF modulator, however, can cause 
the control to interact with the dynamics and result in self-sustained, or limit cycle, 
oscillations. 

A. LIMIT CYCLE DETERMINATION USING DESCRIBING FUNCTION 
ANALYSIS 

In order to predict the existence and characteristics of limit cycles in a system, 
describing function analysis can be employed. Consider the block diagram depiction of a 
system which has been separated into linear and nonlinear elements (Figure 51, adapted 
from Reference 2). 


Nonlinear Element 

Describing Function Linear Element 


r(t) 


^ x(t) 

N(A,(0) 

w(t) 

LiJco) 

y(t) 

J 

i 










Figure 51. Control System with Nonlinear Element 
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For a sinusoidal input to the nonlinear element, i.e. x(t)= Asinfax), the output w(t) of the 
nonlinear element is generally periodic and hence can be represented by a Fourier series 
expansion. The describing function method makes the approximation that only the 
components corresponding to the fundamental harmonic of the expansion, nco with n=l, 
are to be considered. This approximation is based on the assumption that the linear 
element has low-pass properties and hence tends to attenuate the higher harmonics of w(t). 
The describing function of the nonlineai' element is defined as the complex fundamental- 
harmonic gain of a nonlinearity in the presence of a driving sinusoid [Ref. 14]; 


N(A,(o) = 


phasor representation of output component at frequency co 
phasor representation of input component at frequency co 

A 




(26) 


where the magnitude and phase of the describing function are 
\N(A,(0)\ = J^Ja,^ + b,^ 

^ (27) 

ZN(A,co) = 0 = tan"'(aj/fti) 

The describing function N(A, co), is analogous to the frequency response function, 
H(jco), of a linear system. But whereas the linear frequency response function is 
independent of the amplitude of the input signal, the nonlinear "transfer function" or 
describing function depends on both the input amplitude and frequency. In linear control 
loop design, closed-loop poles can be placed to assure a stable response. With a nonlinear 
element, however, the closed-loop pole locations vary with time due to the gain variation 
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with input amplitude. Hence, the root locations in the complex plane can oscillate between 
locations in the left and right-half plane, giving rise to limit cycle behavior. 

By representing the nonlinear element by a describing function approximation, some 
of the same graphical techniques from linear control analysis can be employed to predict the 
limit cycle behavior and to determine stability margins with respect to the limit cycle 
condition. Consider the configuration from Figure 51, assuming that a self-sustained 
oscillation of amplitude A and frequency (O is established and hence the input, r{t) = 0. 
The following relations can be obtained from the figure: 

x = -y 

w = N(A,(o)x 
y = L(j(0)w 


and hence 

y = L(jQ))N(A,co)(-y ), or 
y{l + L(j(o)N(A,co)) = 0 

Since y?K), the limit cycle frequency and amplitude must satisfy the relationship 

A bode plot representation of a system with nonlinearity is depicted in Figure 52. 
The negative inverse of the describing function gain and phase are superimposed with the 
linear plot. The gain and phase margins are redefined with respect to the nonlinear 
boundary, as opposed to the -180 degree phase and 0 dB gain lines. The limit cycle 
frequency and amplitude corresponds to the point which satisfies Equation (28), the point 
of intersection between the nonlinear and linear gain plots. 
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Figure 52. Gain Margin (Gm) and Phase Margin (Pm) for (a) Linear 
System, and (b) System with Nonlinear Element 


B . ANALYSIS OF PWPF MODULATOR 

For the rigid body case, the limit cycle frequency was easily determined by 
simulation. With a flexible spacecraft the dynamic characteristics of the modulator need to 
be considered in order to predict the limit cycle behavior which can result from the 
interaction of the nonlinearity of the PWPF modulator with the control system. The 
dynamic characteristics can be analyzed using the describing function method. 

Although the describing function method is based on the assumption that only the 
fundamental harmonic of the nonlinear response is significant, the residual higher 
harmonics have been considered [Ref. 14]. The corrected-conventional describing function 
uses the root-mean-square value of the first and third harmonics, with the phase of the 
fundamental. A corrected-conventional single input describing function (SIDF) for the 
PWPF modulator is presented in Reference 1, for which the gain and phase of the SIDF 
boundary, N(AmjO)), are determined to be 
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Gain{dB) = 2Q)\og{BIAJ 
Phase{dc^ = -57.3 tan"'(7^ oo) 

where 

h = U„„- U^ff = hysteresis width 

h y(T^(of + \ 

l + expCA/cr^ty))] 


B = 

^3 = 

h = U^„- = hysteresis width 

0 ) = input frequency (rad / s) 

A = minimum pulse width (s) 




( 29 ) 


A block diagram representation of the continuous-time control system is shown in Figure 
53. 


PWPF Modulator Flexible 



Figure 53. PWPF-Modulated Thruster Control Loop 
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Bode plots for the control loop were developed in Matlab using the continuous-time 
dynamics description from Equation (24). The code is included in Appendix C. The 
bode plot for the first five pole-zero pairs of the plant dynamics is shown in Figure 54. 
The gain plot agrees closely with that found in Matrixx using the DBODE command on the 
discrete-time state space model of the plant (see Figure 32). The phase plot from Matrixx, 
however, shows phase lag at the higher frequencies which is not seen in Figure 54. This 
is thought to be another consequence of the discretization process. The possible effect that 
this discrepancy could have is that the predicted phase margin with respect to a limit cycle 
may be smaller than the analysis using the Matlab model would indicate. 



Figure 54. Bode Plot of Flexible Dynamics 
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Figure 55 shows the plant with position and rate feedback control (with K=30) 
superimposed with the SIDF boundary developed using Equations (29). The nonlinear 
boundary intersects the linear plot for the first five modes, with linear control loop gain 
reduction for the sixth system mode frequency sufficient enough that limit cycle interaction 
is not predicted. Hence five possible limit cycle interactions are indicated; only the third 
structural mode limit cycle at = 21 rad/sec was found to be significant, as will be discussed 
shortly. It is interesting to note that in Reference 1 only the first mode pole (representing 
the first system mode of vibration) needed to be considered because the other modes of 
vibration were far removed in frequency from the first and thus far below the SIDF 
boundary. For the flexible spacecraft model in this study, however, the modal density was 
high enough that higher frequency flexible mode limit cycles could be expected. Hence a 
fairly accurate plant model is desired since exclusion of the higher frequency modes in 
initial analysis would fail to predict the higher frequency limit cycles. 



Figure 55. Control Loop Bode Plot with Nonlinear Boundary, K=30. 
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Note that the position and rate feedback effectively provide lead compensation which is 
indicated in Figure 55 by the increase in gain of the higher frequency underdamped pole- 
zero pairs (in addition to the 32 dB contributed by the total gain factor of 42), and by the 
phase lead. The phase lead is desirable as it provides greater phase margin with respect to 
the nonlinear boundary. TTie gain increase however causes the linear boundary to intersect 
with the SIDF boundary at the higher frequencies where limit cycle interaction is of even 
more concern since higher frequency limit cycles are more wasteful of fuel. With the lead 
effect modeled as ts+1, a decrease in the rate feedback gain t would increase the break 
frequency and hence lessen the gain increase seen for the third-through-fifth pole-zero 
pairs. But it follows from Equation 12 (Chapter III) that this would also decrease the 
damping ratio, which may produce undesirable system response (i.e. too much overshoot). 

The stability of the control loop with respect to the third structural mode limit cycle can 
be evaluated by adding pure time delay to the control loop, which introduces phase lag with 
no gain contribution. With a large enough time delay the phase margin is effectively 
eliminated and limit cycling is predicted to occur. Pure time delay in the s-domain is by 
definition; 

G(s) = e-^^ (30) 

where T is the delay time. 

Employing a Taylor series expansion where higher order terms are neglected, an 
approximation to pure time delay with sample and zero-order-hold [Ref. 1] is 

1 - 1 

■■■■ - . = - (31) 

Ts 7’V/12 +75/2-1-1 

Note that for this analysis the use of continuous-time descriptions of the control loop 
elements is founded because the sampling rate used in the discrete-time simulations is 
sufficiently higher than the dynamic response frequencies. 
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Figure 56. Loop Bode Plot with 60 msec Delay, K=30. 
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C. SIMULATION OF LIMIT CYCLE BEHAVIOR 

The analytical stability margin was verified by adding loop time delay to the simulation 
block diagram and performing simulations for various delay times. In real satellite control 
system implementation, the on-board processing capability often limits the sampling rate 
and the speed with which each subsequent control signal command is computed. A multi¬ 
rate sampling scheme was incorporated to simulate the microprocessor sampling interval 
(which contributes to the delay). The block diagram of Figure 57 incorporates the two 
sampling rates and the process delay. The microprocessor sampling interval was chosen to 
be 40 msec. The total time delay is a combination of the microprocessor sampling interval 
and pure computational delay. The PWPF modulator continues to cycle at 100 Hz 
(providing pulse command updates every 10 msec), but receives an input signal every 40 
msec. Hence the modulator sees a constant input for four 10 msec cycles at a time. 


PWPF Modulator 


Process 



Td 


Figure 57. Higher Fidelity Control Loop Model 


The structural mode limit cycle was excited in the simulations first by placing initial 
conditions on the controlled angle {6) and the modal coordinate {qi) corresponding to the 
deformed shape of the satellite and flexible appendage, for whichever mode was closest to 
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the predicted limit cycle frequency. It was found, however, that just about any initial 
condition would trigger the limit cycle in the presence of the proper amount of time delay. 
Figures 58 and 59 shows the time history of the physical and modal coordinates, 
respectively with initial conditions on 6, qj and q 2 , and a 60 msec time delay. A limit 
cycle of frequency 17.3 rad/sec is established after about four seconds, close to the 
predicted frequency of about 20 rad/sec from the analysis of Figure 56. The limit cycle 
amplitude (from the q plot) is of acceptable magnitude, within the pointing accuracy 
"requirement" of +/-. 1 degrees. Note the decay of all modal coordinates except qs, which 
grows in amplitude to the limit cycle condition. The PWPF modulator output shows 
sustained alternating positive and negative thruster firings, characteristic of the limit cycle 
condition. 

The modal coordinate time histories could be converted into physical arm deflections. 
However, the modal magnitudes are somewhat smaller than in the results from Reference 
7, indicating that the physical deflections are not of concern from the standpoint of the 
possibility of structural damage. Also, the main concern here is the control of the central 
body. 
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Figure 58. Initial Condition Response with K=30, 60 msec Delay 



Figure 59. Initial Cond. Modal Response with K=30, 60 msec Delay 
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Figures 60 and 61 show that with a reduced time delay of 20 milliseconds the limit 
cycle is not excited. Note the .12 Hz oscillation in 6 , the frequency of the first cantilever 
mode. It was not possible to excite limit cycles near the fourth and fifth system modes, as 
predicted by the analysis. 



Figure 60. Initial Condition Response, K=30, 20 msec Delay 
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Figure 61. Initial Condition Modal Response, K=30, 20 msec Delay 


The simulation of Figures 62 and 63 was performed with no time delay. Note that 
after the thruster firing stops, the body oscillates at the first system mode frequency of .15 
Hz. An important implication of on-off thruster control can be noted from the modal 
coordinate response. Note that each thruster firing excites the higher modes, and the 
modes ring down between pulses. This points out one of the shortcomings of trying to 
think of an on-off thruster system in linear response terms. Even though the output in 6 
appears to be typical of a sinusoidally forced linear system, each thruster pulse input to the 
plant actually has high frequency content, a.s we know from the trigonometric series 
expansion of a scjuare wave. Note that the sampling rate of 100 Hz is six tunes the highest 
system frequency, so that proper capture of the qe dynamics is assured. 
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Figure 62. Initial Condition Response with K=30, No Delay 
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Figure 63. Initial Condition Modal Response with K=30, No Delay 
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Figures 64 and 65 show that a slew maneuver can still be performed in the presence 
of the third system-mode limit cycle. 



Figure 64. Slew Maneuver with K=30, 60 msec Delay 



Figure 65. Slew Maneuver Modal Response, K=30, 60 msec Delay 
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The analysis was continued by considering a lower PWPF input gain, predicting the 
limit cycle frequency from the bode plot, then adding delay until limit cycling occurred. 
Figure 66 is the linear bode plot with nonlinear boundary for the case K=10. For the 
case with K=30 the fourth and fifth system mode limit cycles were not excitable, hence the 
SIDF boundary is overly conservative. Assuming the same conservatism for the case 
K=10, it was predicted that only the second mode limit cycle would be evident. The 
predicted frequency from the bode plot is 5.5 rad/sec. 



Figure 66. Loop Bode Plot With Nonlinear Boundary, K=10, No Delay 
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With a delay of 280 milliseconds, the gain margin with respect to the second mode 
limit cycle is eliminated (Figure 67). The gain attenuation after 10 rad/sec is due to the 



Figure 67. Loop Bode Plot With Nonlinear Boundary 


K=10, 280 msec Delay 

The predicted stability margin with respect to the second mode limit cycle was verified 
and the results are shown in Figures 68 and 69. The limit cycle frequency is 4.7 rad/sec 
(.75 Hz). Simulations were run with delays of up to 270 milliseconds, but limit cycle 
behavior was not established until 280 milliseconds of delay was added. The assumption 
that the SDDF formulation is conservative was once again verified as the third mode limit 
cycle could not be excited for the case K=10. Note that the limit cycle amplitude violates 
the pointing requirement. 
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Figure 68. Initial Condition Response with K=10, 280 msec Delay 



Figure 69. Initial Cond. Modal Response, K=10, 280 msec Delay 
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This analysis has shown that the control-structure interaction predicted using the SIDF 
formula is overly conservative. The K=30 case predicts limit cycling near the fourth and 
fifth structural modes, and likewise with K=10 a third-mode interaction is still predicted. 
Since the describing function method is by definition an approximate method of predicting 
instability, this is acceptable. But the value of backing up the analysis with simulation is 
clearly demonstrated. 
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VI. CONCLUSIONS 


A pulse-width pulse-frequency modulator has been successfully implemented and 
shown to provide pseudo-linear operation over a certain range of input. Simulations using 
the modulator indicate that this technique provides closer-to-linear actuation than 
conventional bang-bang thruster control methods. This facilitates a more tailorable 
response, which means that a single set of thrusters can be used to meet various control 
objectives. The resulting smoother control results in less thruster firing for identical 
maneuvers as compared to simple bang-bang and time-optimal bang-bang controllers. 

Analysis of the control loop was shown to be important because of the interaction 
between the modulator nonlinearity with the flexible dynamics. Nonlinear describing 
function analysis proved to be useful in predicting flexible-mode limit cycle behavior. 

A. RECOMMENDATIONS FOR FURTHER STUDY 

• The PWPF modulator in this study required a fairly large input gain to provide 
adequate pointing accuracy. The resulting loop gain caused the bandwidth of the controller 
to overlap with flexible modes. The various elements of the modulator could be studied 
and further optimized for the specific system in hand in order to reduce the excitation to the 
structure introduced during maneuvers. Filtering of the first flexible mode (as in Reference 
1) could also be applied. 

• Operation of the PWPF modulator and comparisons to the performance of other 
controllers should be verified by experiments with the Flexible Spacecraft System. 

• Use of thruster pulse modulation for momentum dumping should be investigated with 
the goal of reducing the angular offset noted in Reference 7. 


75 


• A combined-controller approach should be investigated. For instance, the minimum 
time solution could be implemented to provide initial slewing, with switching to a PWPF- 
modulated controller to avoid excessive thruster activity at the termination of the maneuver. 

• The combination of center body control with active control of the flexible appendage 
using piezoelectric actuators and sensors is recommended. The "smart structures" 
experiments to date have been performed with the appendage cantilevered from a fixed 
center body. A multiple-input multiple-output controller could be developed which 
simultaneously considers end point control of the flexible appendage with the center body 
pointing control. 
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APPENDIX A 


The following graphs depict the modulator output for various constant inputs. 
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The following is a listing of the digital implementation of the PWPF modulator. The 
code was implemented within the format of the general Userblock subroutine 
(USROl.FOR). Refer to Figure 5 (Chapter 2) for definition of the parameters. The 
Userblock is implemented in simulations by first compiling the progam ($ FOR USROl) 
followed by linking ($ MWSLINK USROl) which creates an appended MWS executable 
file with default name MYMWS. Entrance into Systembuild is then accomplished using the 
command $MWS MYMWS. See Chapter 5 of the Matrixx/Systembuild Modeling and 
Simulation manual for more information on Usercode Blocks. 


c INITIALIZATION 

IF(TIME.EQ.O)THEN 

c Modulator parameters 
K=4.5 
Tau=. 12976 
lm=.54 
um=l 
uon= .45 
uoff= .15 
dt= .01 

c limiter 

if(al.gt.lm) al= Im 
if(al.lt.-lm) al= -Im 

c +fire (schmidt trigger) 
if(al.gt.uon)then 
ml7= 1. 
else 

ml 7=0. 
end if 

if(al.gt.uoff)then 
nl7=l. 
else 
nl7=0. 
end if 

if((ml7+nl7).gt.O.) then 
ol7=l. 
else 
ol7=0. 
end if 
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c -fire (schmidt trigger) 
if(al.lt.-uon)then 
pl7=l. 
else 
pl7=0. 
end if 

if(al.lt.-uoff)then 
ql7=l. 
else 
ql7=0. 
end if 

if((pl7-i-ql7).gt.0.)then 
rl7=l. 
else 
rl7=0. 
end if 

c initial output 

a3= um*(ol7-rl7) 

y(2)= a3 // y is output vector of Usercode Block 
go to 100 
END IF 

c MAIN LOOP 

c new input u from simulation 
a2= u(l) 

c filter 

al= dt/Tau*K*(a2-a3) + alold*(l.-dt/Tau) 
c limiter 

if(al.gt.lm) al=lm 
if(al.It.-Im) al=-lm 

c + fire (schmidt trigger) 
if(al.gt.uon) then 
ml8=l. 
else 

ml 8=0. 
end if 

if(al.gt.uoff)then 
nl8=l. 
else 
nl8=0. 
end if 

if((ml8-t-nl8*ol7).gt.0.)then 
ol8=l. 
else 
ol8=0. 
end if 

c - fire (schmidt trigger) 
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if(al.lt.-uon)then 
pl8=l. 
else 
pl8=0. 
end if 

if(al.lt.-uoff)then 
ql8=l. 
else 
ql8=0. 
end if 

if((pl8+ql8*rl7).gt.0.)then 
rl8=l. 
else 
rl8=0. 
end if 

c OUTPUT 

a3= um*(ol8-rl8) 

yO)=al 

y(2)=a3 

c reassign the fire commands and filter output for next time step 
ol7=ol8 ^ 

rl7=rl8 
alold=al 

100 RETURN 
end 
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The foUowing listing of a sample simulation/ plotting routine is included for the benefit 
of the next Matrixx user. These routines help minimize the amount of time spent retyping/ 
modifying commands on the Matrix^ command line. The program is accessed using 
oWedt filename, and executed using oexecute ('filename'). In the editor use Ctrl z to 
get back to the editor command line, 'quit' to exit without saving revisions, 'exit' to exit 
and save revisions, 'ch' to get back to editing. 

dt=.01; 

tmax=100; 

t= [0:dt:tmax]'; // Define time vector 

11= 1; // lower limit for graphs 

ul= 6000; // upper limit for graphs 

u= ones(t)*10*pyi80; // input for simulation (10° step) 

y= sim(t, u); // simulation command 

plot(t(ll:ul),y(ll:ui,4),’upper left, ylabeyposition(deg)/ xlabel/time(sec)/') 
plot(t(ll:ul),y(ll:ul,3),'upper right, ylabeyrate(dps)/ xlabeytime(sec)/ date') 
plot(y(:,4),y(:,3),'lower right, xlabel/posittion(deg)/ ylabel/ang vel/... 
title/phase plane trajectory/') 
pause 

hard //hardcopy of plot 
\\print/queue=ln03r matplot.ps 

plot([t(ll:ul) t01:ul) t01:ul)l, [y01:ul,9) y(ll:ul,10) y(ll:ul,l 1)].... //strip chart plot 
'strip ylabel/qllq2lq3/ xlabel/timeCsec)/ ymin-1.2 ymaxl.2') 
pause 
hard 

\\print/queue=ln03r matplot.ps 
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Discrete Super-Block Sanipling Interval First Sample Ext.Inputs Ext.Outputs Enable 

PWPF modulated Ctrl multirate_O.QlQQ 0. 1 12 Parent 


APPENDIX B 
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[simple bang-bang with deadband] 








Discrete Super-Block Sampling Interval First Sample Ext.Inputs Ext.Outputs Enable 

pwpf early implementation 0.0156 0. 1 4 Parent 


Attempts were made (unsuccessfully) to implement the PWPF modulator 
using a Usercode Block for the Schmidt trigger and Systembuild elements 
for the remainder of the loop. 
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APPENDIX C 


These Matlab programs were used for development of the loop transfer function and 
single input describing function (SIDF) boundary for analysis in the frequency domain. 

% The zeros of the plant transfer function coirespond to the 
% eigenvalues of the appendage only, the cantilever frequecies (rad/sec). 
z(l)=.76; z(2)=2.1; z(3)= 16.7; 
z(4)= 22.3; z(5)= 38.6; z(6)= 104.4; 

% Construct the plant transfer function numerator polynomial: 
zeta= .004; % .4 percent modal damping assumed 
numplant= [1 2*zeta*z(l) z(l)'^2]; 
denfactor= z(l)^2; 
fori=2;6 

numplant= conv(numplant,[l 2*zeta*z(i) z(i)^2]); 
denfactor= denfactor*z(i)^2; 
end 

% The poles of the plant transfer function correspond to the free-free or 
% system modes (rad/sec); 
p(l)= .96; p(2)= 2.6; p(3)= 17.2; 
p(4)= 22.8; p(5)= 39; p(6)= 106; 


% Construct the plant transfer function denominator polynomial: 
denplant= [1 2*zeta*p(l) p(l)^2]; 
numfactor= p(l)^2; 
for i=2:6 

denplant= conv(denplant,[l 0 p(i)''2]); 
numfactor= numfactor*p(i)^2; 
end 
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% Factor in the rigid body dynamics and the product-of-poles''2 to 
% product-of-zeros^2 ratio: 

Iest= 11.4; % estimate of the total MOI 
denplant= Iest*conv( denplant,[ conv([l 0],[1 0]) ]); 
factor= numfactor/denfactor; 
numplant= numplant*factor; 

% Poles and zeros of transfer function (no control) 
axis([-10 10-150 150]) 

rlocus(numplant,denplant) %title(' Poles and zeros of plant xfer function') 

pause 

clg 

axis 


% Bode plot of plant 
w= logspace(-l,2,400); 
bode(numplant,denplant,w), pause 


% Now factor in the controller/feedback loop to form the loop xfer function 


K= 10; 

k=K* 1.3985; 
tau= 4.7746; 
numH= [tau 1]; 
denH= 1; 


% PWPF modulator input gain 
% forward loop gain 
% factor for rate gain 
% H is feedback loop 


numloop= k*conv(numplant,numH); % combine plant and controller 
denloop= conv(denplant,denH); 

bode(numloop,denloop,w) % Bode of loop xfer funtion, GH 
pause 

[maglinear,phaselinear,w]= bode(numloop,denloop,w); %linear mag and phase 
%[Gm,Pm,Wcg,Wcp]=margin(maglinear,phaselinear,w); %linear margins 


% Determine the SDDF boundary (magnitude and phase): 
[magsidf,phasesidf]= sidf; % calls function sidf 
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magsidf= magsidf+1.5; % corrected to reflect [Bayloc] experimental results 


% Plot both linear mag and phase and nonlinear boundary: 
subplot(211),semilogx(w,magsidf,'--',w, 20*logl0(maglinear)),grid 
%title('Linear magnitude with SIDF boundary NO DELAY 24feb') 
%gtext('Nonlinear Boundary') %gtext('Gm') 
xlabel('frequency(rad/sec)'),ylabel('magnitude(dB)') 
subplot(212),semiiogx(w,phasesidf,'-',w,phaselinear),grid 
%title('Linear phase with SIDF boundary') 

%gtext('Nonlinear Boundary') 
gtext('Pm') 

xlabel('frequency(rad/sec)'),ylabel('phase(deg)') 

pause 

% The pure time delay causes a phase lag which shrinks the nonlinear phase 
% margin. Use sample and zero-order-hold formula from [Wie & Plescia] 

T= .280; % 60 msec, delay 
numdelay=[l]; 
dendelay= [T^2/12 T/2 1]; 

% augment numloop and denloop and calculate new magnitude and phase: 

numloop=conv(numloop,numdelay); 

denloop=conv(denloop,dendelay); 

[maglinear,phaselinear,w]=bode(numloop,denloop,w); %linear mag and phase 

% Plot both linear mag and phase and nonlinear boundary again for comparison: 
clg 

subplot(21 l),semilogx(w,magsidf,'-',w, 20*logl0(maglinear)),grid 

gtext('Nonlinear Boundary') 

xlabeK'freq (rad/sec)’),ylabel(’magnitude (dB)') 

subplot(212),semilogx(w,phasesidf,'-',w,phaselinear),grid 

gtext('Nonlinear Boundary') 

gtext('Pm ') 

xlabeK'freq (rad/sec)'),ylabel('phase (deg)') 
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SIDF subroutine: 


function [magsidf,phasesidf]= sidf 
% Modulator characteristics and SIDF analysis from [Wie] 

% Static characteristics: 

Uon= .45; 

Uoff=.15; 
h= Uon-Uoff; 

del=.01; %minimum pulse width 
Um= 1; 

Tm=. 12976; 

Kra= 4.5; 

% Calculate SIDF boundaiy gain and phase: 
w=logspace(-1,2,400); 

Am= (Uon-h./( 1 + exp(del/Tm./w)) ).*( (sqrt(Tm.*w) ).^2 + 1 )/Km; 
Bl= 4/pi*Um*sin(del.*w/2); 

B3= 4/3/pi*Um*sin(3*del.*w/2); 

B= sqrt(B1.^2+B3.^2); 

magsidf= -20.*logl0(B./Am); 

phasesidf= -180-( -57.3*atan(Tm.*w)); 
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